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I.  INTRODUCTION 

This  report  documents  a  Fortran  subroutine  called  CHEMEQ  designed 
to  solve  sets  of  ordinary  differential  equations  of  the  form: 


dn . 
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L  n 
i  1 


(1) 


Here  Q  is  the  formation  rate,  L.n.  is  the  loss  rate,  and  n  is  the 
i  ii  i 

density  of  the  ith  species.  Often  the  time  constants  1/L^  for  the 
various  species  differ  by  many  orders  of  magnitude  and  strong  coupling 
between  species  may  be  present.  If  this  is  the  case,  the  set  of  equa¬ 
tions  (1)  is  considered  "stiff"  and  does  not  lend  itself  readily  to 
numerical  solution  by  classical  methods. 

Subroutine  CHEMEQ  was  developed  to  apply  a  specialized  numerical 
technique  "The  Selected  Asymptotic  Integration  Method" (SAIM)  to  this 
class  of  equations.  The  method  has  a  very  low  computational  overhead 
associated  with  it  and  is  particularly  useful  when  combined  with  a 
transport  algorithm  such  as  the  "Flux  Corrected  Transport"2  module  to 
form  reactive  flow  models.  In  such  applications  computer  memory  is  at 
a  premium  because  copies  of  the  chemical  species  variables  are  required 
at  every  grid  point.  Since  CHEMEQ  is  a  single-step  algorithm,  multiple 
copies  of  the  data  from  several  successive  timesteps  need  not  be  saved. 
Further,  since  CHEMEQ  is  simple-  and  single-step,  no  start-up  penalty 
such  as  evaluating  a  large  Jacobian  matrix  is  exacted  at  the  beginning 
of  an  integration  step.  This  is  also  very  important  because  a  reactive 

flow  application  requires  millions  of  chemistry  integration  start-ups. 
Note:  Manuscript  submitted  August  1  1,  1979. 
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Whenever  the  hydrodynamic  processes  in  the  problem  change  the  vari¬ 
ables,  the  chemistry  calculations  must  be  reinitialized. 

The  efficiency  of  CHEMEQ  is  achieved  by  limiting  the  actual  inte¬ 
gration  to  second-order  accuracy  to  minimize  auxiliary  storage  and 
start-up  expense.  In  reactive  flow  applications,  however,  the  reaction 
rates  are  seldom  known  to  better  than  10%  and  the  hydrodynamics  calcu¬ 
lations  are  seldom  accurate  to  better  than  1%.  Thus  integration  of 
the  chemistry  to  better  than  1  part  in  103  or  104  is  an  expensive  folly. 
In  this  regime  CHEMEQ  seems  to  beat  the  classical  methods  by  about  a 
factor  of  50-100  in  speed  on  test  problems  where  start-up  is  not  a 
consideration.  When  a  coupled  hydro  application  on  many  grid  points 
is  attempted  using  a  parallel  processor,  up  to  three  orders  of  magni¬ 
tude  improvement  seems  possible. 

The  SAIM  method  has  been  applied  successfully  to  such  reactive 
flow  problems  as  high  altitude  nuclear  burst  phenomena,3'4  the  solar- 
induced  ionosphere,5  laser-generated  plasma  interactions, 6  and  the 
chemical  kinetics  associated  with  combustion7  problems.  It  has  also 
been  used  successfully  for  chemical  model  development,  particle 
deposition  in  the  ionosphere,8  and  other  physical  problems  where  stiff 
ordinary  differential  equations  arise. 

The  subroutine  CHEMEQ  is  written  in  simple  standard  Fortran  but 
makes  extensive  use  of  the  pipeline  architecture  of  the  NRL  Texas 
Instruments,  Advanced  Scientific  Computer  (ASC) .  The  subroutine  is 
easily  adapted  to  other  machines  without  loss  of  efficiency. 

A  new  subroutine  is  being  prepared  called  VSAIM  (Vectorized  Selec¬ 
ted  Asymptotic  Integration  Method) .  This  subroutine  applies  the 
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asymptotic  method  used  in  CHEMEQ  to  several  independent  sets  of 
equations  (i.e. ,  grid  points)  simultaneously,  and  thereby  it  takes 
full  advantage  of  the  parallel  processing  capability  of  the  ASC.  This 
subroutine  is  particularly  useful  for  solving  the  chemical  kinetics 
associated  with  hydrodynamic  applications  on  computers  that  have 
parallel  processing  capabilities.  VSAIM  will  be  documented  in  a 
subsequent  report. 

Sections  IX  and  III  describe  the  algorithm  and  its  implementation. 

Appendix  A  describes  the  application  of  CHEMEQ  to  various  problems. 

Appendix  B  gives  the  fortran  listing  of  the  subroutine  together  with 
tables  of  internal  and  argument  list  variable  definitions.  Appendix  C 
gives  the  results  of  a  sample  atmospheric  test  problem  using  CHEMEQ 
together  with  program  listings  which  illustrate  the  application  of 
CHEMEQ. 

II.  ALGORITHM 

CHEMEQ  integrates  a  set  of  coupled  ordinary  differential  equations 
(which  may  include  "stiff"  terms)  of  the  form  (1)  by  a  one-step  algo¬ 
rithm.  The  method  has  very  low  overhead  since  all  that  is  required  to 
start  a  new  integration  step  are  the  current  values  of  the  variables 
and  the  derivatives.  A  second  order  predictor-corrector  method,  which 
takes  special  notice  of  those  equations  determined  at  the  beginning  of 
the  step  to  be  stiff  is  employed  to  continue  the  integration  process. 

The  asymptotic  integration  method  applied  to  the  stiff  equations 
best  treats  the  situation  where  the  solution  is  slowly  changing  or 
nearly  asymptotic  yet  the  time  constants  are  prohibitively  small. 

| 
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This  occurs  when  the  formation  rates  and  loss  rates  are  large,  nearly 
equal,  and  there  is  strong  coupling  between  the  equations.  Thus  the 
stiff  equations  are  treated  with  a  very  stable  method  which  damps  out 
the  small  oscillations  caused  by  the  very  small  time  constants.  If, 
however,  the  formation  rates  and  loss  rates  are  small  compared  to  the 
function  size,  the  simple  classical  methods  can  be  utilized  for  these 
equations  to  give  the  combined  method. 

The  predictor-corrector  algorithm  provides  enough  information  to 
choose  the  subsequent  timestep  size  once  convergence  has  been  achieved. 
For  efficiency  an  initial  timestep  is  chosen  which  approximates  the 
timestep  that  will  be  determined  after  convergence  of  the  predictor- 
corrector  scheme.  This  initial  trial  timestep  is  chosen  independently 
of  the  stiffness  criterion  and  is  determined  such  that  none  of  the 
variables  will  change  by  more  than  a  prescribed  amount.  If  the  forma¬ 
tion  rate  is  much  larger  than  the  loss  rate,  it  is  reasonable  to  assume 
that  and  will  remain  relatively  constant  for  large  changes  in  n . . 
Often  the  initial  change  in  n^  may  be  large  enough  to  equilibrate  the 
formation  and  loss  rates.  Thus  the  initial  trial  timestep  is  chosen 
in  two  ways  as  follows; 

<5t  =  tmin[n./n.,  or  (if  Q.  »  L .  n  . )  1/L  .  ]  (2) 

■'ll  1  11  iJ 

Here  C  is  a  scale  factor,  the  same  value  as  the  convergence  criterion 
described  in  Eq.  (6) .  The  minimum  is  taken  over  the  whole  set  of 
equations.  The  timestep  chosen  by  Eq.  (2)  may  be  larger  than  some  or 
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all  of  the  equilibration  times,  in  which  case  the  corresponding  equa¬ 
tions  would  be  classified  as  stiff.  Nevertheless,  when  solved  by  the 
asymptotic  method,  this  timestep  ensures  that  accuracy  can  be  main¬ 
tained.  When  a  stiff  equation  is  close  to  equilibrium,  the  changes  in 
the  functional  values  over  the  timestep  will  be  small  even  though  the 
adjustment  rate  toward  equilibrium  can  be  very  much  shorter  than  the 
timestep.  When  the  stiff  equation  is  far  from  a  dynamic  equilibrium, 
the  timestep  should  be  scaled  down  proportionally  to  the  equilibration 
time  to  ensure  that  the  transition  to  equilibrium  will  be  followed 
accurately.  This  readjustment,  because  of  the  very  fast  rate,  gener¬ 
ally  takes  place  very  rapidly  after  which  much  longer  timesteps  may  be 
taken. 

After  a  timestep  has  been  chosen,  all  of  the  equations  are  sepa¬ 
rated  into  two  classes,  stiff  and  normal,  according  to  the  criterion. 

L . t  <1  Normal 

i 

or  (3) 

L  t  £  1  Stiff 

i 

where  the  value  of  i  is  problem-dependent  and  is  chosen  by  the  user  to 
invoke  asymptotic  treatment  as  necessary.  In  addition,  the  user  may 
force  asymptotic  treatment  on  any  percentage  of  the  set.  Equation  (3) 
is  applied  first.  Then,  beginning  with  the  equations  with  the  shortest 
characteristic  time  (1/L^)  not  already  chosen  by  application  of 
equation  (3) ,  additional  equations  are  selected  with  increasing  time 
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constants  until  the  percentage  of  the  set  specified  is  satisfied.  If 
the  equation  is  considered  stiff  at  the  start  of  the  integration  step, 
it  is  treated  as  stiff  until  the  step  has  been  completed.  The  two 
types  of  equations  are  then  integrated  by  separate  predictor-corrector 
schemes  but  using  a  simple  asymptotic  formula  to  replace  the  usual 
second-order  corrector  equation  for  all  those  equations  which  were 
determined  to  be  stiff. 

The  predictor  part  of  the  step  is  performed  as  follows: 


n. (1)  =  n.(0)  +  5tF . ( 0)  (Normal) 

11  l 

and  (4) 

6tF. (0) 

ni(1)  =  V0)  +IT6t^To)  (stlff) 

i 

where  (0)  =  Fi[t(0),  n^O)  J.  Here  we  start  at  t  =  t(0)  and  wish  to 
find  n^[t(0)  +  iSt]  ^  n.  (1). 

If  we  let  the  integer  in  the  parentheses  denote  the  iteration 
number  then  n^(k)  is  the  kth  iterated  value  of  n.,  or  an  approximation 
to  ni[t(0)  +  St  ].  The  zeroth  iteration,  n^O),  is  the  initial  value 
at  t(0)  and  n^(l)  is  the  result  of  the  predictor  step.  Also  note  that 
F^(k)  =  £\[t(0)+5t,  n^(k)  J  for  the  derivatives.  The  corrector  formulas 
for  the  two  types  of  equations  are: 


n.  (k+1)  =  n. (0)  +  ~  [F . (0)  +  F.  (k)  j  (Normal) 

1  l  2  l  l 


and 


2<5t[Q.  (k)  -L  .  (0)  n  .  (0)+F.  (0)  ] 
ni(k+1)  =  ni(0)  +  -4  At[OOtL1(0)3  1 - 


(5) 
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By  comparing  n^ (k+1)  with  (k)  on  successive  iterations  using 
the  relative  error  criterion  €  to  satisfy  the  following  equation, 


1  i  a  =  max 


[|ni(k+l)  -  n^k)  |  "I 

n.(k+l)6  J  ' 


(6) 


the  convergence  of  each  of  the  individual  equations  can  be  determined. 

As  applied  in  CHEMEQ,  6  is  typically  10“ 3  and  if  the  formation  and 
loss  rates  are  nearly  equal  a  will  be  scaled  down  slightly.  This 
allows  quicker  convergence  for  equations  that  are  nearly  in  equilibrium. 

In  practice  n.  is  constrained  by  a  minimum  value  when  n^  is  decay¬ 
ing  exponentially  toward  zero.  This  lower  bound  is  chosen  by  the  user 
and  must  be  carefully  selected  to  insure  that  its  value  in  no  way 
affects  the  physics  but  yet  decouples  the  equation  from  accurate  inte¬ 
gration.  Decoupling  is  accomplished  by  skipping  the  application  of 
equation  (6)  to  all  equations  that  have  decayed  to  values  correspon¬ 
ding  to  their  lower  bounds.  Convergence  for  these  equations  is  then 
trivial  and  the  function  no  longer  affects  the  size  of  the  timestep. 

For  equations  that  are  decaying  exponentially  to  zero  with  time 
constants  that  are  small  enough  to  control  the  timestep, it  is  important 
for  efficiency  reasons  to  decouple  these  equations  at  the  largest  lower 
bound  possible.  However,  it  must  be  remembered  that  spurious  results 
may  occur  in  other  equations  sensitive  to  the  limited  equations  if 
their  lower  bounds  are  too  large.  This  results  because  the  value  of 
the  function  after  decoupling  is  frozen  at  the  lower  bound  for  the  dura¬ 
tion  of  the  integration  process  or  until  the  total  rate  becomes  posi¬ 
tive.  If  there  is  any  question,  it  is  better  to  be  on  the  conservative 
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side  by  choosing  the  minimum  values  smaller  than  necessary.  This  may 
result  in  a  little  loss  of  efficiency  but  will  reduce  the  possibility 


of  erroneous  results. 

We  have  found  that  maximum  speed  is  realized  by  keeping  the 
allowed  number  of  iterations  on  the  corrector  small.  We  typically  use 
one  or  two.  If  satisfactory  convergence  of  all  equations  has  not  been 
obtained  before  or  during  the  last  iteration,  the  step  is  started  over 
with  a  smaller  timestep.  By  keeping  the  maximum  number  of  iterations 
small,  a  minimum  amount  of  time  is  wasted  on  an  unstable  or  nonconver- 
gent  step  only  to  find  out  that  the  iteration  procedure  did  not  con¬ 
verge.  By  the  same  token,  we  have  found  it  best  to  reduce  the  time- 
step  sharply  (a  factor  of  2  or  3)  when  nonconvergence  is  encountered 
rather  than  to  reduce  it  slowly.  Less  time  is  wasted  this  way  getting 
down  to  a  sufficiently  small  step  for  convergence  if  the  initial  esti¬ 
mated  step  size  is  found  to  be  too  large.  On  the  other  hand,  when 
increasing  the  timestep,  as  for  example  when  convergence  is  achieved 
on  the  first  or  second  iteration,  we  have  found  it  best  to  only  in¬ 
crease  by  5-10%  each  step.  During  the  integration  of  several  succes¬ 
sive  steps,  we  use  the  appropriately  modified  timestep  from  one  con¬ 
verged  integration  cycle  as  the  trial  timestep  for  the  next  integra¬ 
tion  cycle  rather  than  using  Eq.  2.  The  timestep  modification  is  per¬ 
formed  as  follows 


St  =  St  +  .  005^ 


(7) 
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Using  o  as  the  starting  value,  the  /a  is  estimated  with  three 
iterations  of  Newton's  method.  This  gives  the  desired  asymmetrical 
property  in  that  &t  decreases  faster  than  (St  would  increase  for  the 
inverse  value  of  a.  in  addition,  <St  is  modified  very  little  when  o  is 
near  1. 

Once  convergence  of  all  the  equations  is  achieved,  the  new  values 
of  the  n^(6t)  are  set  equal  to  the  values  of  n^(k+l).  One  can  obtain 
convergence  and  completion  of  an  integration  step  after  only  two  deri¬ 
vative-function  evaluations  even  when  some  or  all  of  the  equations  are 
stiff . 

III.  HOW  TO  USE  CHEMEQ 

The  Selected  Asymptotic  Integration  Algorithm,  as  described  in 
Section  II,  has  been  coded  in  Fortran  which  may  be  implemented  on  any 
digital  computer  of  moderate  size.  It  is  intended  as  a  very  fast  but 
moderately  accurate  integrator  which  can  be  used  at  each  grid  point  of 
a  large  hydro-  or  magnetohydrodynamic  calculation.  Single  point  calcu¬ 
lations  are  easily  and  efficiently  accomplished  by  CHEMEQ  as  well. 

CHEMEQ  has  four  entries  which  are  available  to  perform  the  vari¬ 
ous  aspects  of  the  integration.  The  main  entry  is  used  for  normal 
operation.  The  others  provide  flexibility  and  optional  controls.  The 
variables  in  the  argument  lists  and  internal  variables  are  documented 
in  detail  in  Appendix  B. 

1.  CHEMEQ  (DTCHEM,  DFE,  N,  F,  FMIN)  advances  the  equations  the 


specified  interval  DTCHEM. 


2.  CHEMSP  (EPSMN,  EPSMX ,  DTMN ,  TNOT,  PASS,  TASS,  PRT)  resets 
the  specified  control  parameters  if  the  default  values  are  not 
satisfactory. 

3.  CHEMCT  (TMK)  is  for  information  purposes.  This  entry  prints 
information  which  indicates  how  efficiently  the  integration  pro¬ 
cess  has  been  since  either  the  last  call  to  CHEMSP  or  the  last 
call  to  CHEMCT. 

4.  CHEMPR  is  for  diagnostic  purposes.  This  entry  may  be  called 
whenever  an  error  occurs  which  can  be  attributed  to  the  results  of 
CHEMEQ .  A  partial  set  of  the  internal  variables  is  printed  as  a 
diagnostic. 

CHEMEQ  is  the  main  entry  and  is  called  to  advance  the  equations 
as  required.  The  initial  values  are  passed  in  as  arguments.  After 
being  advanced  by  the  integration  they  are  passed  back  in  the  same 
place.  One  of  the  arguments  of  CHEMEQ  is  the  name  of  the  derivative 
function  subroutine  DFE,  utilizing  a  useful  feature  of  Fortran  which 
gives  the  user  the  option  of  specifying  various  configurations  for  the 
derivative  functions  within  the  confines  of  the  same  problem. 

CHEMSP  is  called  whenever  any  or  all  the  default  values  of  the 
control  parameters  in  the  argument  list  are  not  satisfactory.  Vari¬ 
ables  such  as  the  initial  value  of  the  independent  variable,  the 
absolute  minimum  timestep  allowed,  control  parameters  for  convergence 
of  the  predictor-corrector  combination,  and  the  control  parameters 
which  affect  the  use  of  asymptotics  may  be  reset  here. 

CHEMCT  is  called  for  diagnostic  put  poses.  It  displays  information 
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on  the  numbers  of  derivative  function  evaluator  calls  and  the  number 
of  times  asymptotics  were  employed.  It  also  gives  the  number  of  times 
the  integration  step  had  to  be  restarted  with  a  smaller  timestep  due 
to  lack  of  convergence  of  the  predictor-corrector  scheme  since  the 
last  call  to  CHEMCT,  CHEMSP,  or  since  the  beginning.  This  information 
can  be  very  helpful  in  determining  the  relative  efficiency  of  the 
integration  process. 

CHEMPR  is  called  for  diagnostic  purposes.  If  an  error  in  the 
integration  process  is  suspected,  the  user  may  call  this  entry  to 
print  out  some  of  the  internal  variables.  The  current  values  of  the 
production  rates  [C(I)],  loss  rates  [D(I)],  functions  [F(I)J,  inverse 
time  constants  [RTAU  (I)  ],  total  rates  j_CMD],  estimated  timestep  re¬ 
quired,  from  the  previous  step  the  total  rates  [DFS(I)],  the  functional 
values  [FS(I)  Jr  the  initial  functional  values  [FO(I)  j  and  the  minimum 
values  [FmIN(I)  ]  are  printed  for  diagnostic  purposes. 

Two  subroutines  are  referenced  from  CHEMEQ . 

1.  DFE(F,  C,  D,  T)  The  Derivative  Function  Evaluator  which  cal¬ 
culates  the  derivatives  {rb  }  as  required. 

2.  CHEMER  in  the  subroutine  that  is  called  whenever  CHEMEQ  deter¬ 
mines  that  an  error  has  occurred. 

DFE,  the  derivative  function  evaluator,  must  be  supplied  by  the 

user  to  provide  on  request  by  the  integrator,  the  current  derivatives 

{n  }.  It  is  important  to  note  that  nearly  all  of  the  computer  time 
i 

spent  in  the  integration  process  for  most  problems  is  spent  in  this 
user-supplied  routine.  It  is  therefore  >-'xt  rente!  a  valuable  to  put  the 


extra  effort  into  optimizing  this  routine,  especially,  when  it  will 
be  incorporated  into  a  large  hydrodynamic  code.  Here  are  some  sugges¬ 
tions  for  coding  which  may  help  produce  efficient  operation. 

1.  Avoid  all  unnecessary  repetitive  calculations.  Quantities 
which  can  be  calculated  once  should  be  stored  for  subsequent  use. 
In  particular,  divisions  and  mathematical  functions  should  be 
avoided  since  these  are  costly  operations  on  most  machines. 

2.  Replace  all  complicated  functions  with  table  look-ups  when¬ 
ever  possible.  This  can  be  a  very  important  economy  measure. 

3.  Arrange  the  code  in  a  fashion  which  takes  advantage  of  your 
computer's  optimization  features.  For  example,  the  use  of 
register  to  register  or  parallel  processing  capabilities. 

4.  The  user  may  often  take  advantage  of  the  structure  of  the 
problem  he  is  working  with. 

For  example,  in  a  large  atmospheric  reactive-flow  hydrodynamic  code, 
the  density  values  may  vary  drastically  from  the  bottom  of  the  grid 
to  the  top.  Often  at  the  top  a  much  simpler  reaction  scheme  will  be 
sufficient  to  describe  the  chemistry  than  in  the  middle  or  lower  por¬ 
tions  of  the  mesh.  Here  the  user  may  specify  various  configurations 
of  the  reaction  scheme  appropriate  to  the  grid  region  and  save  a  sig¬ 
nificant  amount  of  computation.  There  are  other  ways  to  improve  effi¬ 
ciency  but  they  may  not  be  as  obvious  as  these  listed.  Often  with  a 
little  imagination  and  persistence  combined  with  a  thorough  knowledge 
of  the  problem  area,  significant  improvements  in  computational 
efficiency  can  be  realized. 
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CHEMER  is  called  whenever  CHEMEQ  determines  that  a  severe  error 
has  occurred.  Currently  the  only  error  which  can  be  identified  by 
CHEMEQ  is  when  the  timestep  becomes  too  small.  CHEMEQ  at  this  point 
provides  output  that  may  be  useful  and  then  calls  CHEMER.  The  default 
version  of  CHEMER  does  nothing  but  print  a  message  indicating  that 
CHEMER  has  been  called  and  then  stops  execution.  However,  the  user  may 
supply  his  own  version  of  CHEMER  which  could  provide  printout  of  a 
much  more  complete  set  of  diagnostics  than  CHEMEQ  does  or  manipulate 
the  data  in  such  a  fashion  that  the  integration  process  might  proceed. 

The  actual  arguments  and  internal  variables  used  in  CHEMEQ,  its 
entries  and  the  associated  subroutines  will  be  described  in  detail  in 
the  appendix  sections  of  this  report. 

IV .  SUMMARY 

CHEMEQ  is  intended  to  be  a  general  purpose  integrator  for  a 
specific  type  of  equations.  It  employs  a  very  low-overhead,  moder¬ 
ately  accurate,  low-order  technique.  To  obtain  results  for  most 
physical  models  with  an  acceptable  degree  of  accuracy,  CHEMEQ  can  be 
extremely  efficient.  In  many  areas  where  problems  arc  so  computation¬ 
ally  expensive  they  seem  impossible  to  do  by  other  methods,  CHEMEQ 
gives  accurate  results  in  a  reasonable  amount  of  time.  CHEMEQ  can 
also  be  employed  in  the  development  of  chemical  or  mathematical  models 
when  efficiency  is  not  so  important,  but  the  user  should  not  expect 
eight  figures  of  accuracy.  Two  or  three  figures  over  a  long  integra¬ 
tion  is  a  more  realistic  estimate.  CHEMEQ' s  forte  lies  in  the  solu¬ 
tion  of  the  stiff  ordinary  differential  equations  associated  with 
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chemically  reactive  flow  problems.  Here  the  reaction  rates  are  split 
off  from  the  hydrodynamic  part  of  the  equations  and  solved  separately 
for  each  hydrodynamic  timestep  and  at  each  grid  point.  The  moderate 
accuracy  of  the  methods  used  to  solve  the  hydrodynamic  equations 
suggest  that  the  application  of  a  more  sophisticated  technique,  rather 
than  a  low-order,  low-overhead  method  like  CHEMEQ,  would  waste  valu¬ 
able  computer  time  and  could  possibly  render  the  problem  so  computa¬ 
tionally  inefficient  that  it  would  be  impractical  to  pursue. 

A  potential  user  must  be  aware  that  CHEMEQ  is  not  user-proof, 
problem-independent  and  can  not  always  be  used  as  a  black  box.  The 
method  is  not  identially  conservative  for  arbitrarily  large  timesteps 
when  asymptotics  are  employed  and  the  minimum  values  should  be  chosen 
with  some  thought  since  they  can  become  sources  of  spurious  errors  if 
not  chosen  small  enough  initially. 

All  methods,  such  as  the  selected  asymptotic  integration  method, 
which  do  not  conserve  particle  density  or  charge  balance  automati¬ 
cally  may  be  forced  to  do  so  by  at  least  two  techniques.  In  one  tech¬ 
nique,  conservation  can  be  restored  by  adding  the  various  concentrations 
to  find  the  errors  and  then  by  distributing  these  errors  throughout  the 
densities  in  a  number-conserving  manner.  The  major  fault  with  this  is 
that  a  portion  of  the  errors  is  incorporated  into  concentrations  from 
which  the  errors  may  not  have  arisen.  The  second  and  better  method  is 
to  reduce  the  frequency  of  the  asymptotic  treatment  or  decrease  the 
timestep  size  to  the  point  where  errors  due  to  nonconservation  are 
within  tolerable  limits.  Significant,  improvement  in  computational 
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efficiency  still  results. 

CHEMEQ  is  written  in  standard  Fortran  and  should  be  easily  adapt¬ 
able  to  any  computer  that  accepts  Fortran.  Although  the  present  ver¬ 
sion  is  written  in  a  fashion  that  promotes  vector ization  by  the  ASC, 
no  special  features  of  the  ASC  Fortran  were  incorporated  into  the  code. 

The  storage  requirements  of  CHEMEQ  are  proportional  to  the  maxi¬ 
mum  number  of  equations  for  which  storage  has  been  reserved.  For  a 
maximum  of  25  equations  CHEMEQ  requires  about  2000  words  of  memory 
on  the  ASC. 

Since  CHEMEQ  uses  a  convergence-dependent  algorithm  and  an  adap¬ 
tive  timestep,  the  overall  timing  will  be  strictly  problem-dependent. 

It  will  depend  on  such  things  as  the  coupling  between  and  relaxation 
times  of  the  equations.  As  mentioned  before,  most  of  the  integration 
time  will  be  spent  in  the  derivative  function  evaluations  of  which 
there  are  at  least  two  required  per  CHEMEQ  call.  At  least  50  psec  of 
ASC  CPU  time  are  required  as  integrator  overhead  per  integration  step 
per  equation.  This  does  not  count  the  time  required  to  evaluate  the 
derivatives . 

If  CHEMEQ  is  applied  as  intended,  the  subroutine  can  solve  large 
systems  of  stiff  ordinary  differential  equations  more  efficiently  than 
methods  currently  available.  In  some  cases,  its  efficiency  is  un¬ 
rivaled  . 
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APPENDIX  A 


Table  Al.  Logical  Sequence  of  Calls  for  Chemical  Kinetics 
Without  Transport 


Ca 


culate  initial  conditions,  control  parameters,  etc. 


►Start  loop  over  timesteps: 


CALL  CHEMSP 


Vly' 

CALL  CHEMEQ 


Call  only  to  change  default  values 
of  control  parameters. 


Advance  rate  equations  one  timestep. 


Print  diagnostics  as  needed. 


■End  loop  on  timesteps. 


END  when  loop  over  timesteps  is  complete. 
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Table  A2.  Logical  Sequence  of  Calls  for  Chemical  Kinetics  and 
Transport  Combined 


Ca 


culate  initial  conditions,  grid,  control  parameters,  etc. 


— ^St^lrt  loop  over  hydro  timesteps: 


InYok 


f 


e  transport  algorithm 


tart  loop  over  grid  points: 


CHEMSP 


Call  only  to  change  default  values 
of  control  parameters. 


CHEMEQ 


Advance  rate  equations  at  each  grid 
point  one  hydro  timestep. 


loop  on  grid  points. 


t  diagnostics  as  needed. 


loop  on  hydro  timesteps. 


when  loop  over  hydro  timesteps  is  complete. 
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APPENDIX  B 
LISTING  OF  CHEMEQ 


SUBROUTINE  CHEMEQtDTCHEM,  DFE,  N,  F,  FMIN) 

c 

c 

zo*  ********************************** 

CD 

CD  CHEFEO  (DTCHEM,  DFE,  N,  F,  FMIN) 

CD 

CD  ORIGINATORS!  T.R.  YOUNG  AND  J.P.  BORIS  NRL  1971 

CO 

CD  DESCRIPTION!  CHEMEQ  IS  A  SUBROUTINE  WHICH  SOLVES  A  CLASS  OF 
CO  ORDINARY  DIFFERENTIAL  EQUATIONS  TERMED  STIFF.  THESE  equations 
CO  CANNOT  BE  READILY  SOLVED  BY  THE  STANDARD  CLASSICAL  METHODS  THUS 
CD  THE  SELECTED  ASYMPTOTIC  INTEGRATION  METHOD  IS  EMPLOYED  BY  CHEM£C. 
CO  THE  EQUATIONS  ARE  OIVICEC  INTO  TWO  CATAGORIES  BASED  ON 

CD  EQUILIBRATION  TIMES  AND  ARE  INTEGRATED  BY  EITHER  A  LOW  ORDER 

CD  CLASSICAL  METHOD  FOR  THE  EQUATIONS  WHICH  HAVE  LONG  EQUILIBRATION 
CO  TIMES  OR  A  VERY  STABLE  STEP-CENTERED  METHOD  WHICH  HELPS 
CO  PRESERVE  THE  ASYMPTOTIC  NATURE  OF  THE  SOLUTIONS  WHEN 

CD  EQUILIBRATION  TIMES  ARE  VERY  SHORT,  an  ADAPTIVE  STEPSIZE  IS 

CD  CHOSEN  TO  GIVE  ACCURATE  RESULTS  FOR  THE  FASTEST  CHANGING  CUANTITY. 


CO 

THE  ROUTINE  ASSUMES  THAT 

ALL  OF  THE  INTEGRATED  CU*NTITES  AND  THE 

CD 

CD 

CO 

TIME  STEP  ARE 

POSITIVE. 

ARGUMENT  LIST 

DEFINITION! 

CO 

DTCHEM 

REALM 

THE  INTERVAL  OF  INTEGRATION  OR  THE 

I 

CD 

RANGE  OF  THE  INDEPENDENT  VARIABLE. 

CO 

0.0  «■  T  <■  DTCHEM. 

CD 

DFE 

REAL*« 

THE  NAME  OF  THE  DERIVITIVE  FUNCTION 

I 

CD 

evaluator  subroutine. 

CD 

N 

INTEGER 

THE  NUMBER  OF  EQUATIONS  TO  BE 

I 

CO 

INTEGRATED.  AN  ERROR  EXISTS  IF  N  IS 

CO 

GREATER  THAN  ND  SET  BY  THE  PARAMETER 

CO 

statement. 

CD 

F  (N  ) 

REAL*d 

THE  INITIAL  VALUES  AT  CALL  TIME  I/O 

CD 

AND  THE  FINAL  VALUES  AT  RETURN  TIME. 

CD 

FMIN(NJ 

REAL* A 

MINIMUM  VALUES  FOR  EACH  FUNCTION. 

I 

CO 

CO  LANGIAGE  AND  LIMITATIONS!  ALTHOUGH  THIS  SUBROUTINE  IS  WRITTEN  IN 
CO  A  FASHION  WHICH  PROMOTES  VECTORIZATION  BY  THE  ASC  COMPILER ,  THE 
CD  FORTRAN  IS  NEARLY  STANDARD  AND  SHOULD  WORK  WITH  MINOR  MOCIPlCAT- 
CD  IONS  ON  ANY  MACHINE. 

CD 

CO  ENTRY  POINTS!  FOUR  ENTRY  POINTS  ARE  PROVIDED  FOR  FLEXIBILITY  AND 
CD  OPTIMUM  CONTROL. 

CO 

CD  CHEMEQ l  ADVANCES  THE  EQUATIONS  THE  GIVEN  INCREMENT  ' DTCHEM ' , 

CD 
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CD 

CHEf'CT  t 

INFORMATIVE,  PRINTS  THE  VALLES  OP  THE  INDICATIVE 

CD 

COUNTERS  LISTED  9EL0WJ 

CD 

1.  THE  NLFBER  OF  TI“ES  ASYHPTOTICS  HERE  USED. 

CO 

2.  the  NUMBER  of  DERIVATIVE  FUNCTION  EVALUATIONS. 

CO 

3.  t>e  number  TINES  THE  INTEGRATION  step  was  restarted 

cc 

DUE  TO  NONCONVERGENCE  OF  THE  PREDICTOR-CORRECTOR 

CO 

SCHEME. 

CO 

CD 

CHE^SPj 

PROVIDES  THE  USER  WITH  THE  OPTION  TO  RESET  THE  HOST 

CD 

IMPORTANT  CONTPOL  PARAMETERS. 

CD 

CO 

CFEMPRi 

INFORMATIVE,  PRINTS  OUT  INTERNAL  VARIABLES  FOR  DIAGNOSTIC 

CD 

PURPOSES. 

CO 

CO 

SUBROUTINES  REFERENCED! 

CC 

CD 

CFEl 

WHOSE  ACTUAL  NAME  AND  DEFINITION  ARE  SUPPLIED  BY  THE  USER 

CO 

IS  CALLED  to  OBTAIN  THE  DERJVITIVE  FUNCTIONS. 

CC 

CD 

CALL  DFE (F »  C,  C,  T) 

CO 

APGLMEN  T 

LIST  TO  CFEl 

CC 

FCN) 

REAL*"  CURRENT  VALLES  OF  THE  DEPENDENT 

I 

CO 

variable. 

cc 

C  (N ) 

REAL*"  CALCULATED  FORMATION  RATES. 

0 

CO 

CCN) 

REAL*-  CALCULATED  loss  rates. 

0 

cc 

T 

REAL*"  CURRENT  VALLE  OF  THE  INDEPENDENT 

I 

CD 

VARIABLE. 

CO 

CD 

CHEFERi 

IS  CALLED  WHENEVER  AN  ERROR  IS  DETECTED.  CURRENTLY  THE 

CD 

ONLY  ERROR  RECOGNIZED  IS  A  TIME  STEP  THAT  IS  TOO  SMALL. 

CC 

CD 

CALL  CHEMER 

CO 

arglment 

LIST  TO  CHEMERl  NO  ARGUMENTS. 

CO 

CO* 

k 

c 


c 

c 


c 

c 


PARAMETER 

NC  ■  15 

0 

RFAL*e 

TS.  TN 

REAL 

F (N/ND/ 

),  FMIN (N/N  C/  )  , 

CCNTJ, 

D  (ND) 

REAL 

FS(ND). 

DFS(NC),  FC(ND) 

,  SCR  A  £ 

ND), 

REAL 

ASY (ND) 

,  COR (ND ) 

INTEGER 

FCOUNT, 

ACOLNT,  RCOUNT, 

TFCNT , 

TACN 

CATA 

FCOUNT, 

ACOUNT,  RCOUNT, 

TFCNT , 

TACN 

,  RTAUS(NC) 
SCRB(ND),  RT  AL (K  D ) 

T,  TRCNT 
T,  TRCNT/6*0/ 


i 
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DATA  PASYI/O.OC/,  TCRASY/lOO.n/,  NDD/ND/ ,  EPSCL/100.0/ 

DATA  TFC/Z<m000C8/,  DTHIN/1  .OE-15/,  SQREPS/0.50/ 

DATA  TSTART,  DT,  DTTEST/3*0  .CX,  TN/c.0D«0C/,  C/NC*0.0/ 

CATA  EFSMAX/10.0/,  LO/fc/,  EPSM  IN/  1 .0E-O2/ ,  C/ND*0.0/ 

C 

C  TEMPERARY  FIX  I  SEE  DO  LOOP  13C  1  l  1  1 
CATA  SIGNM/ZSOOOOOOO/ 

c 

C  CHECK  INPUT  PARAMETERS. 

IF  CN  .IE.  N  CD )  GO  T8  110 
KRITECLO,  10023  N,  NCD 

1002  F8RMAT(5C/),  i  FROM  -CHEMEQ-  I  NO.  OF  EQ.S  REQUESTED  IS  TOC, 
'  LARGE'/'  REQUESTED  (',15,'),  MAX.  ALLOWED  (',15,')') 

STOP 

C 

c  INITIALIZE  the  control  parameters. 

110  TN  «  O.OC+CO 

CTTARG  ■  OTCHEW 
C 

C  STORE  AND  LIMIT  TO  'FMIN'  THE  INITIAL  VALUES. 

CO  1  I  *  1,  N 
C  ( 1 3  ■  0.0 
cm  ■  o.o 

FOCI)  ■  F  ( 1 3 

1  F ( I )  ■  AWAXICF(I),  FMIN(I)) 

C 

C  EVALUATE  THE  CERIVITIVES  OF  THE  INITIAL  VALLES. 

CALL  OFE(F,  C,  0,  SNGLCTN  +  TST  ART ) ) 

FCOUNT  •  FCOUNT  ♦  1 

ESTIMATE  THE  INITIAL  STEPS1ZE. 

STRONGLY  INCREASING  FUNCTI0N5CC  >>>  C  ASSUMED  HERE)  USE  A  STEP- 
SIZE  ESTIMATE  PROPORTIONAL  TO  THE  STEP  NEEDED  FOR  THE  FUNCTION  TO 
REACH  EQUILIBRIUM  WHERE  AS  FUNCTIONS  DECREASING  OR  IN  EQUILIBRIUM 
USE  A  STEPSIZE  ESTIMATE  CIRECTLT  PROPORTIONAL  TO  THE  CHARACTER¬ 
ISTIC  STEPSIZE  OF  THE  FUNCTION,  CONVERGENCE  OF  THE  INTEGRATION 
SCHEME  IS  LIKELY  SINCE  THE  SMALLEST  ESTIMATE  IS  CHOSEN  FflR  THE 
INITIAL  STEPSIZE. 


SCRTCH 

B 

l.OE-RO 

CO  15  I 

1,  N 

SCRA(I) 

. 1 tEPSM IN  * AB5 (C  ( I ) )  - 

D  ( I ) 

SCRB(I) 

SIGNC1.0/F  Cl), 

SCRA  ( I 

)) 

SCR  A  ( I ) 

SCR8 (I ) *C  C I ) 

scfiecn 

-AQS(ABS(C(I)) 

-  CCI) 

)«SCRBC1) 

SCR  A  ( I ) 

AMAX1 (SCR  A  ( I ) , 

SCR B  C I 

)) 

SCRTCH 

B 

AMAXl (SCRA  (I) , 

SCRTCH) 
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CT  *  SQREP5/SCRTCH 

IF £CT  .IT.  TCRASY)CT  ■  SQRT (DTaTCR ASY ) 

IF  (DT  .ST .  CTCHEHJCT  ■  DTCMEH 
C 

C  THE  STARTING  VALUES  ARE  STORED. 

100  TS  ■  TS 

ASSIGN  ASYHPT8TIC  OR  NORHAL  TREATMENT  FOR  EACH  EOLATION  AT  THE 
BEGINING  OF  EACH  STEP. 

NASY  ■  ACOLNT 

EQUATIONS  WITH  TOO  SHORT  A  CHARACTERISTIC  STEPSIZE  ARE  SELECTED 
FOR  ASYMPTOTIC  TREATMENT. 

CO  130  I  ■  1,N 
RTAU(I)  «  C(I)/F(I) 

FS (I )  «  F  Cl 3 
C F S  (13  •  C(I)  -  0 Cl 3 
SCR  A  (I )  ■  RTAUm  -  TCRASY 

THE  FOLLOWING  TWO  CARDS  REPLACE  THE  THIRD  WHICH  COES  NOT  COMPILE 
PROPERLY  ON  NX  -  5.027.131. 

SCRBfJ)  *  ANDISIGNH,  SCRA(I)) 

ASY(I)  b  ,5  t  OR  ( .5 #  SCRB(I)) 

ASY (1 5  ■  .5  ♦  SIGN ( ,5 #  SCRA(D) 

COR C 1 3  ■  OF S C 1 3  »  6CI3*A9V(S) 

RTAUSCI)  »  PTAU  < I ) * A5Y  (1) 

130  ACOUNT  ■  ACOUNT  ♦  ASY(I) 

NASY  ■  PASYI*N  -  ACOUNT  ♦  NASY 
IF  (NASY  .LT.  1 )G0  TO  101 

COMPLETE  THE  SELECTION  OF  EQUATIONS  FOR  ASYMPTOTIC  TREATMENT  UP  TO 
THE  PERCENTAGE  • P ASY I  •  .  EQUATIONS  WITH  THE  SHORTEST  CHARACTERISTIC 
5TEPSIZE8  ARE  CHOSEN  FIRST. 

CO  20  I  «  1 .NASY 
RTHX  «  -1.0Et7Q 
CO  35  J  «  1  * N 

IF  (ASY (I)  .CT.  O.DGO  TO  35 
IF(RTMX  .GT.  RTALCJnCO  TO  35 
RTHX  •  RTAKJJ 
JS  •  J 

35  CONTINUE 

ASY  C JS )  >  1.0 
RTAUS(JS)  B  RTAU(JS) 

COR  ( JS)  i  COR  ( JS)  -  C  C JS 3 
20  ACOUNT  ■  ACOUNT  ♦  1 

c 

C  FIND  THE  PRECICTOR  TERMS. 
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on  n  oo  ooo  ooo 


m  :o  r  i  s  i,n 

c 

C  THE  FIRST  ORCER  PREDICTION  F5B  the  ASVHFnTIC  FUNCTIONS  EE  Cl  EE  3  T» 

C  ELLER'S  METHOD  e;r  the  N0NASYmET0T:C  FUNCTIONS  Ir  "RUl"  s  0 . " . 

SCEPCI)  a  CFSm/Cl.y  ♦  DT.PTALSCI)) 

5  CONTINUE 

c 

C  L I  h*  I T  crcPEASIt'G  FLNCTI"f>S  TO  THEIR  HTM'ij*  VALUES. 

CO  105  I  s  1  ,N 

SCEPCI)  *  FS(I)  t  CT*SCRGCI) 

105  F  C I )  =  AHAAS  (SCRBCI),  F  H1 1 N  ( I )  J 

TK  s  TS  a  CT 
C 

C  EVALUATE  THE  DtRIVIUVCS  FOR  THE  CORRECTOR. 

CAlL  PFccF,  C#  D,  Sf-GLCTt.  t  TSTArtj) 

FCOIN'.T  3  rc“UNT  +  1 
EPS  =  1.CF-1J 

ce  y  i  s  i ,  r, 

STEP  CENTERED  CORRECTOR  FOR  THE  ASYMPT0TIC  FUjN  C  T I  r« -s  PEDUCCS  TO 
THE  NOTIFIED  EULCR  METHOD  FOR  the  (,0|)ayS,,PT',TIC  functions, 

SCRR(I)  s  ASY  1 1 ) *D  c  I J 
RTALCI)  5  SCRB(I)/F(I) 

y  scppct)  s  ccopci)  ♦  cm  ♦  cscrpcd  -  pern) 

.  /  (2.o  ♦  .s*ct*cptauci)  *  RTAusmn 

CALCUUATC  f-EH  F,  CHECK  fjr  CONVERGENCE,  AnC  LImIT  DECREASING 
PL NOTICES.  THE  ORDER  **F  THE  BPEERATISNS  P'  This  LOOP  IS  IMPORTANT. 
CO  6  I  S  1  ,  N 

scrp(l)  3  A^AVICFSCI)  ♦  DT«sCPeCI),  O.C) 

SCRACI)  S  ASSCSCPRCI)  -  FCI)) 

F ( I )  a  AHAX1  (SCRCCI),  F  M  I  N  ( I  )  ) 

scra(I)  s  scRAcn/Fcn 

SCALE  RELATIVE  ERROR  DO*r.  *hen  C  &  n  ARE  NEARLY  CCUAL. 

SCRPCI)  s  APS(CCI)  -  Clin/CCC)  ♦  r  ci )  ♦  l.OF-30) 

SCRP(I)  S  AHINlCSCfiBCI),  SCRACI)) 

SCRACI)  s  SCPACI)  ♦  SCRB(I) 

REMOVE  cCLATlVt  EpROR  C On TR IUIjT I  ON  IE  FUNCTION  valjf  IS  LESS  Than 
THE  f  IM"U."  vALLE. 

SCR°  C I )  =  ,25*  (FS  ( I )  ♦  F  ( I ) )  -  Ft-  in  c  I ) 

SCRP(I)  s  ,25  ♦  SIGN  ( .25,  SCRACI)) 

SCRA(I)  3  SCHfi  C I ) *SCR  A  ( I ) 
fc  EPS  s  AHAXJ C5CRACI),  EPS) 

EPS  s  ErS*tp5CL 
C 


* 
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C  PRINT  OUT  CIANCSTICS  IF  STEPS  I ZE  BECOMES  TOO  SHALL. 

IF  (DT  .CT.  CTH IN  ♦  1.0E-16*TN)GP  TO  40 
HRITECLO,  1003)  CT.  TN,  DTPIN 
CO  25  L  I  l.N 
CPC  •  CCD  -  0(D 

CTC  a  EP£PIN*F{D/(ABS(CPO)  ♦  1.0E-30) 

25  HRITECLO,  1000)  CCD.  CCD,  F  CD ,  RTAU(L) ,  CPC,  CTC,  CFSCL), 

.  FS  CD  ,  FCCD,  FMIN(L) 

1003  FORHATC'l  CHEPEC  ERROR!  STEPSIZE  TOO  SPALL  1  1  J',  /, 

1  '  DT  ■  1PE1C.3,  '  TN  »  025.15, 

2  '  CTHIN  ■  '  ,E  1 0 .3 ,  //,  1RX,  'C,  9X,  *C,  9X,  *F',  6X,  ’RTAL  1 , 

3  5X,  'C  -  D  CTC  CFS 1 ,  8X,  'FS',  8X,  'FO  F^IN') 

1000  F0RPATC5X.  1P12E10.3) 

CT  «  OTOHEP  -  TS 

CT  a  APINllCTHIN,  ABS  CCT ) ) 

C 

C  CALL  ERROR  DIAGNOSTIC  ROLTINE 

CALL  CHEPER 
C 

C  CHECH  FOR  CONVERGENCE. 

00  IF (EPS  .CT.  EPSHAX)GO  TO  30 

C 

C  END  CHECK. 

CTTARG  ■  DTCHEH  •  TN*TFD 
IF(DTTARG  ,GT.  0.0)GO  TO  1C 
RETURN 
C 

C  PERFCRP  STEPSIZE  MODIFICATIONS. 

30  RCOUNT  a  RCOUNT  ♦  1 
TN  a  TS 
C 

C  ESTIPATE  SORT  CEPS)  BY  NEKTON  ITERATION. 

10  RTEPS  »  .5*  CEPS  ♦  1.0) 

CO  50  J  i  1,2 

50  RTEPS  ■  .5*  CRTEP5  ♦  EPS/RTEPS) 

CT  a  DT*  C 1 ,0/RTEPS  ♦  .005  ) 

CT  a  APIMCCT,  SNGL  CTFC  *  (DTCHEP  -  TN))) 

C 

C  BEGIN  NEW  STEP  IF  PREVIOUS  STEP  CONVERGED. 

IF (EPS  .CT.  EPSPAX)GO  TO  IC1 

CALL  DFECF,  C,  D,  SNGL (TN  +  TSTART) ) 

FCOUNT  a  FCOUNT  ♦  1 
GO  TO  100 
C 
C 

ENTRY  CHEPCT  CTPK) 

C 
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CD* 

t  t  t 

CO 

CO 

CHEFCT 

(THK) 

CO 

WRITE  OUT  THE 

VALUES  OF  THE 

VARIOUS  INDICATIVE  COUNTERS  THAT  THE 

CO 

PROGRAH 

KEEPS 

• 

CO 

CD 

ARGUHENT  LIST 

DEFINITION! 

CO 

TMK 

REAL*4 

A  FLOATING  POINT  NUHBER  PRINTED 

CO 

to  identify  the  call. 

CO 

CD 

CUTPLT 

VARIABLE  DEFINITION! 

CO 

THK 

REAL*4 

FLOATING  POINT  IDENTIFIER. 

CD 

FCOUNT 

INTEGER 

NUHBER  OF  DERIVATIVE  SUBROUTINE  CALLS 

CO 

SINCE  THE  LAST  CALL. 

CO 

ACOUNT 

integer 

NUHBER  OF  TIHES  THE  ASYNPTOTIC  TREAT- 

CD 

MEAT  HAS  USED  SINCE  THE  LAST  CALL. 

CO 

RCOUNT 

INTEGER 

NUHBER  OF  TIHES  STEPSIZE  WAS  REDUCED 

CD 

SINCE  LAST  CALL. 

CD 

TFCNT 

INTEGER 

TOTAL  OF  FCOUNT  TO  THIS  CALL. 

CO 

TACNT 

INTEGER 

TOTAL  OF  ACOUNT  TO  THIS  CALL. 

CO 

TRCNT 

INTEGER 

TOTAL  OF  RCOUNT  TO  THIS  CALL. 

CO 

CO* 

C 

*  *  * 

TFCNT  ■  TFCNT  ♦  FCOUNT 
T*CNT  «  TACNT  *  ACOUNT 
TRCNT  a  TRCNT  +  RCOUNT 
C 

C  PRIM  BUT  INDICATIVE  COUNTERS. 

NR ITE (LO •  1000)  TM«,  FCOUNT,  ACOUNT,  RCOUNT,  TFCNT,  T  ACM , 
.  TRCNT 

1000  FORNATC  CHEMEQ  INDICES)  THK  a  ',  IPE1Q.3, 

.  *  FCOUNT ,  ACOUNT,  RCOUNT  a  ',  317,  '  TOTALS!  ',  317) 

C 

C  RESET  COUNTERS. 

FCOUNT  ■  0 
ACOUNT  ■  0 
RCOUNT  a  0 
RETURN 
C 
C 

ENTRY  CHEMSPIEPSMN,  EP5WX ,  OT*N,  TNOT ,  PASY,  TASY,  PRT) 

C 

CO 
CO 
CO 
CO 


CHEHSPCEPSHN,  EPSHX,  DTHN  ,  TNOT,  PAST,  TASY,  PRT) 


CD 

RESET 

ANY 

LOCAL 

CONTROL  PARAMETERS  IF  THEIR  RESPECTIVE  INPUT 

CO 

VALDES 

ARE 

GREATER  Than  ZERO.  CEFAULT  VALUES  ARE  USED  IF  THE 

CO 

CO 

CO 

INPIT 

VALUES 

ARE 

ZERO  lR 

LESS  REPECTJVELV. 

ARGUMENT  LIST 

DEFINITION! 

CO 

EPSPN 

REAL*« 

THE  MAXIMUM  RELATIVE  ERROR  ALLOWED 

CO 

FOR  CONVERGENCE  OF  THE  CORRECTOR  STEP 

CO 

CEFAULT  VALLE  1  l.QE-02 

CD 

EP3PX 

REAL*« 

THIS  NtHgER  provices  the  basis  for 

CD 

DECIDING  WEATHER  CONVERGENCE  CAN  eE 

CD 

ACHIEVED  WITH  OLT  ADDED  ST£PSIZE 

CO 

REDUCTION,  IF  EPS/EPSHIN  IS  GREATER 

CO 

THAN  EPSMX  FURTHER  REDUCTION  IS 

CO 

APPLIED. 

CO 

CEFAULT  VALLE  I  10.0 

CD 

DTPS 

REAL*« 

THE  SMALLEST  STEPSI2E  ALLOWED. 

CO 

CEFAULT  VALUE!  l.DE-15 

CO 

TNOT 

RE AL*A 

THE  INITIAL  VALLE  OF  THE  INDEPENDENT 

CO 

variable  t. 

CO 

DEFAULT  VALUE!  0,0 

CD 

PAST 

RE  AL*  a 

the  PERCENTAGE  of  the  EOLATIONS  FOR 

CO 

WHICH  ASYMPTOTICS  WILL  ACHATS  8E 

CO 

APPLIED.  EOLATIONS  WITH  THE  SMALLEST 

CO 

CHARACTERISTIC  S7EPSIZE  ARE  CHOSEN 

CO 

FIRST. 

CO 

TASV 

REAL*a 

ASYMPTOTICS  ARE  APPLIED  IF  THE  CHAR- 

CO 

AC7ERISTIC  STEPSIZE  OF  AN  EQUATION  IS 

CO 

LESS  THAN  TASV, 

CD 

DEFAULT  VALUE!  UOE-02. 

CO 

CO 

CO 

CO 

CO* 

PRT 

REAL*R 

controls  the  output  of  chemsp,  any 

NON  ZERO  VALUE  SUPPRESSES  ALL  PRINT 
OUTPUT  FROM  THIS  ENTRY, 

*  *  *  * 

*  * 

ft  ft 

« 

ft  ft  *  ft  ft 

ft*ftftft«ft****ft«ft*4ftftft* 

c 

EPSMIN  »  1  .OE-02 

IF  {EP3HN  ,GT .  Q.OJEPSMIN  ■  EP5MN 
IF (EP3MN  .GT.  O.OTSQREPS  •  S ,0*SGRT (EPSM IN ) 
EPSCL  *  I.C/EPSMIN 
EPSMAX  •  lC.O 

IF(EP5MX  .GT,  C.OEPSMAX  ■  EPSMX 

ITMIN  »  l.CE-15 

IF(CTPN  . GT ,  0,0)DTM1N  *  OTMN 

73T ART  a  T|\0T 

PAST!  •  C.C 

IF (PAST  .GT.  0.O)PASYl  «  ,C1«(PASV  ♦  .5) 

TCP  AST  •  1C0.0 
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IF  (TASY  .GT.  O.OHCRASY  •  J.O/T*SV 
C 

C  PRIM  SEW  VALUES  flF  CONTROL  PARAMETERS. 

IFfPRT  .EG.  0.0) 

.  WRITECL8,  1001)  EPSMN,  EPSMX,  CTpN,  TNOT,  PAST,  TASY 
1001  FORMATt*  IMTAU2E  "ChEPEG"  VIA  -CHEMSP"',  /, 

.  •  EPSMN,  EPSMX,  CTHN,  TNOT,  PASY,  T*SY  ■  ',  lPoG10.3) 

RETURN 

C 

ENTRY  CHEMPR 
C 

Ct!*  *********************************  * 

CO 

CO  CHEMPR  MAY  Pi  CALLED  WHEN  EVER  AN  ERROR  OCCURS  ThAT  CAN  EE 

CD  ATTRACTED  TO  TME  RESULTS  OF  CHEMEG.  A  PARTIAL  SET  OF  THE  INTERNAL 

CD  VARIABLES  IS  PRINTED  AS  A  DIAGNOSTIC. 

CD 

CD*********************************** 

C 

WRITECLO,  1003)  D7 ,  TN,  CTTEST 
CO  05  L  «  l.N 
CMC  •  C  CL  3  -  0  (L ) 

CTC  *  EPSMIN*FCL)/(ABS{CMD)  ♦  l.CE-30) 

05  NRITE(LO,  1000)  C  (L  )  *  CCD,  F  CL ) »  RTAU(L),  CMC,  DTC,  CFSCUi 
FSCL),  F CCD,  FMIN(L) 

c 

RETURN 

END 
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LISTING  OF  CHEMER. 


subroutine  ci-efer 

DIAGNOSTIC  ROUTINE  FOR  STIFP  fl.C.E.  SOLVE*  -CHEMEQ- 
FRIN'T  10C1 

»  F0RHAT{5(/),  '  LIBRARY  VERSION  OF  -CHEMER-  CALLED.',  /, 

•  LSERS  MAY  SUPPLY  TMEIR  OKN  VERSION  FOR  DIAGNOSTICS.',  /, 
'  NO  ARGLMENTS  ARE  REQUIRED.',  /, 

.  '  PROGRAM  MILL  CONTINUE  RESETTING  THE  STEP  SIZE  TO  MIN-', 

•IMUHS  IF  A  NORMAL  RETURN  IS  MADE.',  //, 

.  »  (STOP  69)  EXECUTED  FROM  LIBRARY  VERSION  OF  -CHEMER-') 

STOP  69 

end 
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These  subroutines  may  be  punched  onto  cards  directly  from 
the  program  listing.  On  computers  other  than  the  TI  ASC  the  PARAMETER 
statement  sould  be  removed  and  occurrences  of  ND  in  the  declarations 
should  be  replaced  by  a  fixed  point  number  at  least  as  large  as  the 
largest  set  of  equations  to  be  integrated.  This  subroutine  should  be 
complied  on  the  K  level  optimization  on  the  ASC  unless  the  number  of 
equations  expected  is  small.  Then  the  J  level  compilation  will  give 
the  most  efficient  code.  No  other  compilation  options  are  required 
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Fortran 

Mathematical 

Variable 

Type/Origin 

Variable 

Comments 

ACOUNT 

I/L 

Index 

Counter;  records  the  number 
times  asymptotics  were 
employed. 

ASY(I) 

R/L 

Logical 

Records  the  location  of  the 
equations  selected  for  asymp¬ 
totic  treatment. 

C(I) 

R/L 

Qi 

Current  formation  rates. 

CMD 

R/L 

Q  -L  n 
l  i  l 

Current  total  rate  (inter¬ 
mediate  variable  for  printing) 

COR(I) 

R./L 

Multiple  Usage 

Temporary  storage  array. 

D  (X) 

R/L 

L  n. 
l  l 

Current  loss  rates. 

DFS(I) 

R/L 

Q.-L  n. 

Ill 

Total  rate  saved  from  the 
beginning  of  the  step. 

DT 

R/L 

St 

Current  timestep. 

DTC 

R/L 

fit 

Timestep  suitable  for  stabi¬ 
lity  estimate  (intermediate 
variable  for  printing) 

DTCHEM 

R/A 

0<ts  t  . 

chem 

Range  of  the  independent 
variable  t. 

DTMIN 

R/L 

fit  . 
mm 

Minimum  timestep  allowed 
Default  value;  1.0E-15) 

DTMN 

R/A 

fit  . 
mm 

Minimum  timestep.  Replaces 
DTMIN  if  DTMN  >  0. 

DTTARG 

R/L 

St  -  t 
chem 

Intermediate  variable  used 
for  end  check. 

DTTEST 

R/L 

St  . 

min 

Intermediate  variable  used  to 

check  DT  for  minimum  value. 

EPS 

R/L 

Max  ( o . ) 
i 

The  maximum  value  of  the  rela¬ 
tive  error.  Used  to  check 
for  convergence. 

EPSCL 

R/L 

i/e  . 

mm 

Intermediate  variable  used  to 
avoid  repeated  divisions. 
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Comments 


Fortran  Mathematical 

Variable  Type/Origin  Variable 


EPSMAX 

R/L 

e  /t  . 
max  mm 

If  EPS  is  larger  than  this 
value  the  step  is  restarted. 
(Default  value;  10.) 

EPSMIN 

R/L 

t  . 
mm 

The  convergence  criterion. 

DT  for  following  step  will  be 
scaled  proportional  to  EPSMIN/ 
EPS.  (Default  value;  .01) 

EPSMN 

R/A 

£  . 
mm 

The  convergence  criterion. 
Replaces  EPSMIN  if  EPSMN  >  0. 

EPSMX 

R/A 

t  /c 

max  min 

Step  restart  criterion. 
Replaces  EPSMAX  if  EPSMX  >  0. 

F  ( I) 

R/A 

n . 

1 

The  current  values  of  the 
solution  to  the  set  of  equa¬ 
tions  being  integrated. 

F  0{  I) 

R/L 

n.  (0) 

1 

Initial  values  at  t  . 

o 

FCOUNT 

I/L 

Index 

Counter ,  records  the  number 
derivative  function  calls. 

FMIN(l) 

R/A 

ni (min) 

Minimum  values  for  each 
equation. 

FS(I) 

R/L 

n.  (0) 

i 

The  values  of  the  solution 
saved  from  the  beginning  of 
the  current  step. 

I 

I/L 

Index 

Subscript  counter. 

J 

I/L 

Index 

DO  loop  subscript . 

JS 

J/L 

Index 

Save  location  for  specific  J. 

L 

I/L 

Index 

DO  loop  subscript. 

I/L 


LO 


Numerical  value  for  the 
logical  unit  for  the  printed 
output . 


Fortran 

Variable 


Comments 


Mathematical 


N 

I/A 

The  number  of  equations  to  be 
advanced. 

NASY 

I/L 

Intermediate  used  in  the  asymp¬ 
totic  selection  process. 

ND 

I/L 

Constant 

Array  size  specification  set 
by  the  PARAMETER  Statement. 

This  is  an  ASC  Fortran  feature. 

NDD 

I/L 

Storage  location  for  ND.  This 
is  ASC  specific  Fortran. 

PASY 

R/A 

0sPASY<100 

Percentage  of  equations  to  be 
treated  by  asymptotics.  Re¬ 
place  PASYI  if  PASY  >  0. 

PASYI 

R/L 

% 

The  percentage  value  of  the 
set  of  equations  that  will 
always  be  selected  for  asymp¬ 
totics.  (Default  value;  0.) 

PRT 

R/A 

Print  control  parameter.  If 
non- zero  printer  output  from 
entry  CHEMSP  sill  be  suppres¬ 
sed. 

RCOUNT 

I/L 

Index 

Counter,  records  the  number 
of  times  integration  process 
had  to  be  restarted  due  to 
non-convergence . 

RTAU(I) 

R/L 

L. 

i 

The  reciprocals  of  the  charac¬ 
teristic  times. 

RTAUS(X) 

R/L 

L. 

1 

The  reciprocals  of  the  charac¬ 
teristic  times  saved  from  the 
beginning  of  the  current  step. 

RTEPS 

R/L 

/Max  (a  . )  /£  . 

i  min 

Used  to  estimate  new  time- 
steps. 

RTMX 

R/L 

L. 

i  max 

Intermediate  variable  used  to 
store  the  maximum  value  from 

RTAUS . 

SCRA(I) 

R/L 

Multiple 

usage 

Temporary  storage  array. 

SCRB(I) 

R/L 

Multiple 

usage 

Temporary  storage  array. 
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m 


Fortran 

Variable 

Type/Origin 

Mathematical 

Variable 

Comments 

SCREPS 

R/L 

A  . 

mm 

Square  Root  of  EPSMIN. 

SCRTCH 

R/L 

Intermediate  variable  used  to 
estimate  the  initial  timestep. 

TACNT 

I/L 

Index 

Counter  total.  Records  the 
total  each  time  ACOUNT  is  set 

to  zero. 

TASY 

R/A 

T 

Asymptotic  treatment  selection 
criterion.  Replaces  TCRASY 
if  TASY  >  0. 

TCRASy 

R/L 

T 

Asymptotic  treatment  selection 

criterion.  This  parameter  is 
problem-dependent  and  the  value 
should  be  proportional  to  the 
overall  characteristic  time- 
step  of  the  system  of  equations 
being  solved.  (Default  value; 
0.1.  This  value  is  often  suit¬ 
able  for  high  altitude  atmo¬ 
spherical  and  many  combustion 
problems.)  Often  it  is  useful 
to  vary  this  parameter  as  the 
solution  progresses. 


TFCNT 

I/L 

Index 

Counter  total.  Records  the 
total  each  time  FCOUNT  is  set 

to  zero. 

TFD 

R/L 

Round-off  parameter.  Should 
have  a  5  in  the  last  signifi¬ 
cant  figure  for  single  preci¬ 
sion  floating  point  words. 

TN 

D/L 

t 

Current  value  of  the  indepen¬ 
dent  variable  t. 

TNOT 

R/A 

t 

o 

Initial  value  of  the  indepen¬ 
dent  variable  t.  Replaces 
TSTART  if  TNOT  >  0. 

TRCNT 

I/L 

Counter  total.  Records  the 
total  each  time  RCOUNT  is  set 

to  zero. 
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Fortran 

Variable 

Type/Origin 

Mathematical 

Variable 

Comments 

TS 

D/L 

t 

o 

The  value  of  the  independent 
variable  t  saved  from  the 
beginning  of  the  current  step. 

TSTART 

R/L 

t 

o 

Initial  value  of  the  indepen¬ 
dent  variable  t.  (Default 

value;  0) 

TMK 

R/A 

Floating  point  number  (typically 
the  value  of  Time)  printed  to 
identify  the  call  to  CHEMCT . 

APPENDIX  C 

This  example  involves  the  integration  of  seven  rate  equations 
which  describe  the  time  evolution  of  an  atmospheric  chemical  relaxa¬ 
tion  test  problem  with  cesium  and  cesium  ions.  This  particular  set  of 
rate  equations  which  was  originally  suggested  by  D.  Edelson  of  Bell 
Laboratories  is  considered  stiff  and  not  well  suited  for  numerical 
integration  by  classical  methods. 

The  sample  program  listed  in  this  section  is  designed  to  deter¬ 
mine  the  efficiency  of  various  stiff  ordinary  differential  equation 
solvers  on  this  test  problem.  In  this  example  CHEMEQ  is  used.  Effi¬ 
ciency  is  determined  by  comparing  the  results  at  the  end  of  the  inte¬ 
gration  interval  with  known  values  and  the  computer  time  required  to 
obtain  these  results  or  various  values  of  the  convergence  parameter 


Table  Cl.  A  List  of  the  Seven  Species  Together  with  Their 

Initial  and  Accepted  Final  Concentrations  for  the 
Test  Problem 


i 

Species 

Number  Densities 

Initial 
y±,  cm-3 

Final 

Yi-  cm-3 

i 

e~ 

1.0 

X 

102 

4.9657897283 

X 

104 

2 

°2~ 

5.2 

X 

102 

2.5913949444 

X 

104 

3 

Cs+ 

6.2 

X 

102 

7.5571846728 

X 

104 

4 

Cs 

1.0 

X 

1012 

1.5319405460 

X 

103 

5 

Cs02 

0 

1.000  X  1012 

6 

N2 

1.4 

X 

101  5 

1.400  X  101:; 

7 

°2 

3.6 

X 

1014 

3.590  X  1014 

Table 

C2.  A  List  of  the  Seven  Reactions  and  Reaction  Rates 
Through  which  the  Seven  Species  of  the  Test 
Problem  Interact 

No. 

Reaction 

Rate  constant  or 
frequency 

1 

02~  +  Cs+  ■*  Cs  +  02 

5 

X  10" 8  cm3  s_1 

2 

Cs+  +  e~  -*■  Cs  +  hv 

1 

X  10-J2  cm3  s"3 

3 

Cs  +  hv  -*■  Cs+  +  e~ 

3. 

24  X  10-3  S"1 

4 

02“  +  hv  -*  02  +  e" 

4 

X  10-1  s"] 

5a 

02  +  Cs  +  M  ->  Cs02  + 

M 

1 

X  10“ 3  *  cm6  s-1 

6 

02  +  e-  +  02  -*•  02-  + 

°2 

1. 

24  x  lO-38  cm6  S-1 

7 

aM  = 

02  +  e“  +  N2  -*■  02"  + 
[Cs]  +  [Cs02]+  [N2j+  [0. 

N2 

2  J> 

1 

X  10“ 31  cm6  s_1 

In  the  following  listing  TACSR  is  the  main  program  which  provides 
the  logic  and  overall  control.  Initialization  and  output  of  results 
takes  place  here.  CSDFE  is  the  derivative  function  evaluator  for  the 
text  problem.  Results  for  nine  values  of  the  convergence  parameter 
EPS  are  printed  at  the  end  of  the  section. 
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PROGRAM  T*C3F 

THIS  IS  THE  EXECUTIVE  PROGRAM  THAT  PROVIDES  THE  LOGIC  NESSICARY 
TO  ACVAKCE  A  REACTIVE  SEVEN  SPECIES  TEST  PROBLEM  POR  AN  EVALUATION 
OP  THE  INTEGRATION  METHOD  FOR  VARIOUS  VALUES  OP  THE  CONVERGENCE 
PARAMETER.  IN  THIS  EXAMPLE  "CHEMEQ*  HILL  BE  EMPLOYED. 

PROGRAM  SPECIFICATIONS. 


REAL*0 

DSEC 

REAL 

V  !  1 0 ) ,  YP(10),  YMINUO),  Yl  C 10) ,  EPSILI10), 

EPS !  15) 

INTEGER 

SPSYMI10) 

external 

CSCFE 

DATA 

OATA 

DATA 

YMIN/lOAl.OE-OO/,  MXCASE/9/ 

SP8YM/  '  02-  ' ,  •  C  S  ♦  ' ,  'CS',  '  CS02  ' ,  '02',  'N2' 
EPS/.l,  .05,  .01,  .005,  .001,  .0005,  .0001, 

,  'NE'/ 
.00005, 

.00001/ 


1000  F0RMATC1CASE  NO.  15,  *  PARAMETERS!',  /, 

.  •  CONVERGENCE  PARAMETER  EPS  »  ',  1PE10.3,  /, 

.  •  INNER  LOOP  LENGTH! 15) 

1001  FORMAT!/,  '  SPECIE  Y  -  INITAL  Y  -  PINAL 

.  '  Y  -  SOLUTION  REL  ERR') 

1002  PORMATCSA,  AO,  1P3E15.6,  E10.3) 

1003  FORMAT!/,  •  T  -  INITIAL  •  (',  1PE10.3,  ')  T  -  FINAL  »  !', 

.  E10.3,  ')') 

1000  FORMAT!/'  INTEGRATION  STATISTICS!') 

1005  FORMAT!'  CPU  TIME  USED  FOR  INTEGRATION!',  1PE10.3, 

.  •  SEC.,  CFU  TIME  NORMALIZED! ',  E10.3) 

1006  FORMAT!'  SUM  OF  THE  RELATIVE  ERRORS  SQUARED!  ',  1PE10.3) 

1007  FORMAT!/) 


INITIALIZE  CONTROL  PARAMETERS. 

"TSCALE"  IS  A  NORMILIZATION  FACTOR  USED  TO  COMPARE  EFFICIENCY  OF 
INTEGRATION  CODES  FROM  DIFFERENT  COMPUTER  INST ALATIONS.  "TSCALE" 
MAY  EE  DETERMINED  BY  TIMING  A  TEST  COOE  ON  ALL  INSTALATIONS 
INVOLVED. 

TSCALE  ■  1.0 

C  SET  INNER  LOOP  LENGTH.  SEE  COMMENTS  eELOH  FOR  DEFINITION. 

INLP  ■  1 
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c 

c 

c 


c 

c 


c 

c 

c 

c 


c 

c 


c 

c 


c 

c 


c 

c 


c 

c 


c 

c 


c 

c 


SET  THE  TOTAL  NUMBER  OF  SPECIES  "NS"  AND  THE  NUMBER  TO  BE 

integrated  "na«. 

NS  a  7 
M  ■  5 

"TI»  -  INITIAL  TINE#  "TP*  -  FINAL  TINE. 

T I  ■  0.0 

TF  a  100C.C 

CELTAT  a  tTF  -  TD/INLP 

STORE  INITIAL (T I  ■  0.0)  AND  FINAL (TF  a  1000.0)  VALUES. 


02- 

YI(1) 

» 

S.iOOE+Oe 

YF(1) 

• 

2.59139492C6lDt04 

CS* 

YIC2) 

a 

6.t00E+02 

YF  (2) 

a 

7  #557 18 46 03000 tOO 

cs 

YIC3) 

a 

l.COOE+12 

YF  (3) 

a 

1.53l9«o51722D*03 

C302 

TI{«) 

a 

1.000E-30 

YI(«) 

a 

l.COQE+04 

YF(«) 

a 

9. 999999235 16D*  11 

02 

YI  (5) 

a 

3.600E+14 

YF(5) 

a 

•3. 5900000005 ID  +  19 

N2 

YI  (0) 

a 

1.400E+15 

YF  (0 ) 

a 

1.400000000000*15 

NE 

YI<7) 

s 

1.000E+02 

YF  (7) 

• 

4.965789082390*09 

LOOP 

1  OVER 

THE  TEST  CASES. 

CO  30 

ICASE  a  l.NXCASE 

PRINT 

10CO,  ICASE,  EPS(ICASE), 

CALL  CHENSP  (EPS(ICASE),  0.,  0.,  TI,  0..  10. 0,  0.) 
CPLT  a  o.O 
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c 

C  RESET  •¥■  Tfl  INITIAL  VALLES  "VI». 

C6  35  I  •»  l.NS 
35  Yd)  ■  YId) 

C 

C  SET  TIMER. 

CALL  SECONCd.  DSEC  ) 

INNER  LOOP  TO  DETERMINE  OVERHEAD  OR  RELATIVE  STARTING  EPFECIENCV 

OF  itegration  scheme  being  tested. 

CO  5  ISTEP  ■  t.INLP 
C 

C  CALL  INTEGRATOR. 

CALL  CHEMEC  (DELTAT ,  CSDFE,  NA,  Y,  YMN) 

5  CONTINUE 

C 

C  CALCULATE  CPU  TIME  USED  IN  THE  INTEGRATION  PROCESS . 

CALL  SECOND  (Q»  DSEC) 

CPLT  a  CFUT  ♦  DSEC 
TNORH  a  CPLT/TSCALE 
C 

C  RESET  ELECTRON  DENSITY. 

Y  (7 )  a  VC2)  -  YCI) 

C 

C  CALCLLATE  RELATIVE  ERROR. 

CO  10  I  a  1 .NS 

10  EPSIL(I)  ■  ABSCYCI3  •  YF(I))/AMIM  (Y(I)  ,  YF(I)) 

SUM  a  O.C 
CO  25  I  ■  l.NS 

25  SUM  a  SUM  ♦  EPSIL(I)«*2 
C 

C  PRIM  RESULTS. 

PRINT  10C3,  TI.  TF 
PRINT  1001 
CO  15  2  -a  l.NS 

15  PRINT  1002,  SPSYMCI),  YId),  YFd),  Yd),  EPSIL(I) 

PRINT  10CO 
PRINT  1006,  SUM 
PRINT  1005,  CPUT,  TNORM 
PRINT  1007 
CALL  CHEMCTCTF) 

30  CONTINUE 

STOP  60 

END 
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SUBROUTINE  CSCFECY,  C,  0,  7) 


C 

cc 

CP 

CO 

CD 

CO 

CD 

CO 

CD 

CO 

CO 

CO 

CO 

CD 

CO 

CO 

CO 

CO 

CD 

CO 

CO 

CO 

CO 

CO 

c 

c 

c 


c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


CSOFE  C  Y  #  C,  C.  T) 

DESCRIPTION! 

DERIVATIVE  FUNCTION  EVALUATOR  (CFE >  FOR  AN  ATMOSPHERIC  CHEMICAL 
RELAXATION  TEST  PROBLEM  INVOLVING  CESIUM  ANO  CESIUM  IONS.  FORMAT¬ 
ION  AND  LOSS  RATES  ARE  CALCULATED  FOR  THIS  SET  OF  "STIFF  ORDINARY 
DIFFERENTIAL  EQUATIONS"  THAT  NAS  SUGGESTED  BY  BY  0.  EDEL50N  OF 
BELL  LABORATORIES. 

ARGUMENT  LIST  DEFINITIONS! 

Y  ( 1 1  R«  <1  CURRENT  VALUES  OF  THE  FUNCTIONS  PLUS  THE  I/O 

EXTRA  DATA  at  THE  END  OF  THE  ARRAY  THAT  MAY  BE 
PASSED  BACK  AND  FORTH  BETWEEN  "C3DFE"  AND  THE 
MAIN  PROGRAM.  LOCATIONS  IN  Y(I)  WHICH  REPRESENT 
THE  FUNCTIONS  BEING  ADVANCED  SHOULD  NOT  BE 
TAMPEREO  WITH  HERE. 

CCD  R*R  TOTAL  FORHATION  RATES.  I 

CCD  R*a  total  LOSS  RATES.  I 

T  R*<l  THE  VALUE  OF  THE  INDEPENDENT  VARIABLE.  I 


LOCAL  SPECIFICATIONS. 


REAL  NE,  NJ 

REAL  VCD#  CCD,  0(1) 

UTILIZE  LOCAL  STORAGE  FOR  VARIELES. 

02M  *  V(U 

C3P  ■  Y  1 2 ) 

CS  ■  YC3) 

C 9 0 2  ■  Y(«) 

02  ■  Y (5) 

N2  •  Y (6  J 

CALCULATE  ELECTRON  DENSITY  FOR  LOCAL  USE  ANO  TRANSMISSION  BACK  TO 
THE  MAIN  PROGRAM  VIA  Y(7 J .  HOWEVER  IN  THIS  CASE  THIS  VALLE  SHOULD 
NOT  eC  TRUSTED  SINCE  "CHEHGE*  WILL  NOT  CALL  THE  "DFE "  WITH  THE 
LATEST  FUNCTION  VALUES  AFTER  THE  FINAL  STEP  HAS  CONVERGED.  V (7) 
WILL  BE  ONE  ITERATION  BEHIND  IN  THIS  CASE.  YD)  AND  YCfe)  ARE 
EXAMPLES  THO,  OF  HOW  DATA  MAY  EE  TRANSFERED  BETWEEN  THE  "CFE"  AND 
THE  MAIN  PROGRAM, 

NE  ■  AMAXUCSP  -  02H,  0.0) 
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Y  (7 )  ■  KE 
C 

C  CALCLLATE  REACTION  RATES. 

CR1  ■  5.00£-08*02»*CSP 
CR2  »  t,COE- 12»CSP*NE 
CP3  a  3.e«E-03*CS 
CP4  a  4. COE-Ol *02M 

CRS  a  l.C0E*31«02*CS«(C8  4  CS02  ♦  N2  ♦  02) 

CPt  a  1.24E-3O*02«O2*NE 
CP7  a  1. COE-31 *02*N2*NE 

CALCLLATE  TOTAL  FORMATION  RATES  (CCD)  AND  TOTAL  LOSS  RATES  (D  (D ) 


c 

FOR 

EACH 

SPECIES 

• 

t 

c 

02M 

CM) 

• 

CR6 

♦ 

CR7 

C(l) 

■ 

CR1 

♦ 

CR4 

c 

c 

CSt 

C  (2) 

■ 

CR3 

C(2) 

9 

CR1 

♦ 

CR2 

c 

c 

CS 

CC3) 

m 

CR1 

♦ 

CR2 

C  (3 ) 

9 

CR3 

4 

CRS 

c 

c 

CS02 

C  (4) 

9 

CR5 

c 

c 

02 

C  (5 ) 

9 

CR1 

4 

CR4 

C(5) 

9 

CF5 

4 

CR6  4  CRT 

C 

RETURN 

END 
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m  NO.  1  M««Ki1EI<3l 
OkVERCENCE  MMKCTM  EPS  ■  1  .OOCE'Ol 

NNE*  IOOF  lENGTHI  1 
NITHI2F  •C‘‘E»'E8«  VI*  "CFEl-SP* 

PSHK,  EFSHX,  CTVN,  TKBT .  P*»V,  T*8V  ■  0.100  0.000  0.000  0.000  0.000  10.00 
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CHEMEO  INCICEtl  TPK  •  l.OOOE  OJ  FCSINT,  *£OU*T,  RCOUNT  »  *2*  J05  0  WM.SI  *5J  «T0 


#»*  f  ‘OS  3St3 


iNTECRiTiet*  srmsttcsi 


m  NB.  T  P»R»HETERSI 
0NVER6EHCE  PrRAHETER  EPS  *  I.OOCE-04 
NN£»  l DOF  LESSTHI  I 


INTE6RPTJ0N  STPTISTICSl 

SUP  BP  THE  REL*TIVE  ERRORS  SSUARECl  2.R*4E>03 

CPU  TtHE  ISEC  FOR  IHT£CR*T10H|  6. 12SE  BO  Sec.i  CPU  TIH£  KBRHAL12EDI  6.120E  00 


